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ABSTRACT 

I describe DESPOTIC, a code to Derive the Energetics and SPectra of Optically Thick 
Interstellar Clouds. DESPOTIC represents such clouds using a one-zone model, and can 
calculate line luminosities, line cooling rates, and in restricted cases line profiles using 
an escape probability formalism. It also includes approximate treatments of the other 
dominant heating and cooling processes for the cold interstellar medium, including 
cosmic ray and X-ray heating, grain photoelectric heating, heating of the dust by 
infrared and ultraviolet radiation, thermal cooling of the dust, and collisional energy 
exchange between dust and gas. Based on these heating and cooling rates, DESPOTIC 
can calculate clouds' equilibrium gas and dust temperatures, and their time-dependent 
thermal evolution. The software is intended to allow rapid and interactive calculation 
of clouds' characteristic temperatures, identification of their dominant heating and 
cooling mechanisms, and prediction of their observable spectra across a wide range 
of interstellar environments. DESPOTIC is implemented as a Python package, and is 
released under the GNU General Public License. 
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1 INTRODUCTION 

With the advent of powerful radio telescopes such as the 
Atacama Large Millimeter Array (ALMA), it has become 
possible to study the cold interstellar medium (ISM) in un- 
precedented detail and at greater distances than ever before. 
Observations from these facilities have stimulated a great 
deal of theoretical interest in the properties of the cold ISM, 
both nearby and in environments far-removed from those 
found near the Sun. One of the goals of these theoretical in- 
vestigations has been to study how the thermodynamics of 
gas, and thus the nature of the star formation process within 
it, varies as a function of environment. A second goal has 
been to predict the observable emission of gas in a variety 
of environments. 

Theoretical investigations of this sort often benefit from 
approximate calculations using idealized geometries that can 
produce relatively fast results, while also including a wide 
range of microphysical processes in order to determine which 
ones are important. Exampl es of this sort of code include 
cloudy JFerland et al.l ll99ot ) and dusty jlvezic fc ElitzuJ 



Il997f ). which are widely used in the ISM community. How- 
ever, there are far fewer publicly-available tools capable of 
performing similar functions for the dense, optically thick 
phase of the interstellar mediu m. Traditional p hotodisso- 
ciation region (PDR ) codes (e.g. iMeiierink fc Sp aans 2005; 
iLe Petit et al.1120061 ) perform calculations in this regime, but 
the complexity of the problem means that these codes are 
too computationally-costly for either broad surveys or quick, 
interactive scans of parameter space. A number of authors 
have released codes capable of performing fast calculations 
of molecular emission line spectra using the large velocity 
gradient or various other forms of the escape probability 
appro ximat ion (e.g. CASSI^j and RADEX, Ivan der Tak et al.l 
[2007]). While these are useful tools for the analysis of ob- 
servations, they are only capable of predicting line emission 
given fixed physical conditions, and they do not calculate 
many quantities of interest for theoretical modeling, such 
as rates of heating and cooling, thermal equilibria, or time- 
dependent thermal behavior. 

The need for codes that are capable of performing 
calculations of this sort is apparent from the wide vari- 
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ety of applications they have f ound in th e recent liter- 
ature. For example, I Goldsmith! l|200ll ) and iLesaffre et al.l 
(2005) investig ate the temperature str u cture within pro- 
tostellar cores. I Krumholz fc Thompsonl (|2007l ) use an es- 
cape probability model to study the relationship between 
star forma tion rates and emission is a vari ety of molec- 
ular lines. iKrumholz. Lerov fc McKeei (|2011h use thermal 
equilibrium models of the ISM to explore the relation- 
ship betweeji_star_Jormation and thech«micaAstate_of the 
gas. INaravanan et all J201ll. I2012TI. [Shettv et al.l (|2011al lbh, 
and iFeldmann. Gnedin fc Kravtsovl I 2012al lbl) all investi- 
gate the conversion between observed CO luminosity and 
molecular mass using simulations of galaxies, coupled to 
post-processing to predict the observable line emission. 
INaravanan fc Dave! |2012a lbi) perform calculations of inter- 
stellar medium temperatures as a way of estimating the 
Jeans mass in molecular clouds, and its possible implica- 
tions for changes in t he ste llar IMF over cosmologi cal times. 
iPapadopoulosT |2010l ) and iMeiierink et all |201ll ) consider 
star formation in extreme environments with X-ray and 
cosmic-ray fluxes far higher than are found in the Solar 
neighborhood, and in the process rely on calculations of the 
thermal behavior of gas un der these conditions. Similarly, 
iMunoz fc Furlanettol J2013h study the ISM in high-redshift 
galaxies where the metallicity is much lower and the cosmic 
microwave background is much hotter than in the present- 
day universe. 

With a few exceptions, all of these authors developed 
their own custom codes to model the thermodynamics and 
line emission of the cold ISM. However, this effort is largely 
duplicative, since these calculations all involve the same re- 
lated set of problems: given a set of physical conditions 
(e.g. density, extinction, velocity dispersion) in some region 
of the ISM, calculate the observable line emission, identify 
the dominant heating and cooling processes, and calculate 
either the gas cooling rate or the equilibrium temperature. 
Moreover, the results of the calculations can be difficult to 
compare due to the differing assumptions and approxima- 
tions made by the various authors in their modeling, not all 
of which are well-documented in the literature. 

In order to support theoretical investigations facing 
problems of this sort, reduce duplication of effort, and 
encourage calculations with documented, open-source 
tools to allow easy comparison between authors, I have 
developed a software library to Derive the Energetics and 
SPectra of Optically Thick Interstellar Clouds (DESPOTIC). 
DESPOTIC uses an escape probability formalism to calculate 
line emission, and couples this to a calculation of either 
equilibrium or time-dependent gas and dust temperatures, 
including the dominant processes in a wide variety of 
environments: cosmic-ray and X-ray ionization heating, 
photoelectric heating, grain-gas energy exchange, and 
radiative heating and cooling of dust grains. The software 
is implemented as Python package, enabling easy, inter- 
active calculation, and also easy integration with other 
software. It also provides an automated interface with 
the Leiden Atomi c and Molecular Database (LAMDA; 
Schoier et"ail [2005) . DESPOTIC is publicly available from 



http: //www.ucolick. org/~krumholz/codes/despotic/ , 



https : //code . google . com/p/despotic/source/checkout 



and the Python Package Index, and is released under the 



GNU General Public License. It comes with extensive 
documentation, including a ~ 40 page User's Guide. 

In the remainder of this paper, I describe the model 
system that DESPOTIC uses and the equations it can solve 
(§ [2} and the numerical methods by which it solves those 
equations (§|3J. I then provide some example applications 
(§ IH, provide some warnings about the limitations of the 
code (§[5j), and summarize (§|6}. 



2 MODEL SYSTEM AND EQUATIONS 
SOLVED 

2.1 Physical Model 

The basic physical system treated by DESPOTIC is a uni- 
form spherical cloud (though other simple geometries are 
provided as options, as described below). Such a cloud is 
characterized by several physical and chemical properties, 
which are taken to be uniform unless stated otherwise. The 
physical properties are a volume density of hydrogen nu- 
clei riH and a mean column density of hydrogen nuclei Nu, 
the gas temperature T g , the dust temperature Td, the non- 
thermal velocity dispersion <jnt, and (optionally) a bulk ra- 
dial velocity gradient dv r /dr. Note that DESPOTIC defines 
the column density as an average over the cloud, i.e. it is 
the total number of hydrogen atoms in the cloud divided by 
the cloud's cross-sectional area. 

The dust within a cloud is characterized by six quan- 
tities. Three of these describe the dust cross-section per H 
nucleus to thermal radiation at temperature T = 10 K, <7d,io, 
to radiation in the range of 8 — 13.6 eV that dominates pho- 
toelectron production, <Jd,PE, and averaged over the diffuse 
interstellar radiation field (ISRF) o"d,iSRF. The fourth quan- 
tity is the total dust abundance normalized to the Milky 
Way value, Z' d . The remaining two quantities are the dust 
spectral index j3 for thermal radiation, and the gas-grain col- 
lisional coupling coefficient qgd ■ I define all of these terms 
in detail below. 

DESPOTIC parameterizes the radiation field (including 
cosmic rays) around the cloud by the following quantities: 
C, gives the primary ionization rate per H nucleus due to 
hard x-ray photons and cosmic rays, x describes the energy 
density, normalized to the Solar neighborhood value, of the 
non-thermal interstellar radiation field produced primarily 
by stars, T ra d,dust gives the infrared radiation field seen by 
the dust, and Tcmb is the cosmic microwave background 
temperature. 

Finally, the chemical composition of the cloud is de- 
scribed by the abundances of bulk constituents and trace 
emitting species. The abundances of the bulk constituents 
in the DESPOTIC model are given by xm, i p h 2 , £oH 2 > ihc, 
x c , and x H +, which describe atomic hydrogen, para-Hb, 
ortho-Hb, helium, free electrons, and free protons, respec- 
tively^ All abundances are measured by number density 
relative to the number density of hydrogen nuclei; note that 



2 Although DESPOTIC includes free protons and electrons, it is 
only intended for use in regions where the gas is predominantly 
neutral, i.e. £ H + <C xm + 2(x p n 2 +a; H,), and similarly for x e . 
It does not include many heating and cooling processes that are 
important in highly ionized regions. 
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this means that fully molecular composition corresponds to 
£ p h 2 + x h 2 — 1/2 rather than 1. The abundances of emit- 
ting species (e.g. CO, HCN, H2O, etc.) are characterized in 
the same way, with Xi representing the abundance of the ith 
emitting species. 

Given the bulk composition, one can also compute a 
number of additional quantities, of which we will make use 
below. Three of these are the mean mass per H nucleus /xh, 
the mean mass per free particle /1, and the isothermal sound 
speed c s , given by 



Mh = xm + x u + + 2(x p h 2 + x u 2 ) + iXUc 
M = 



(1) 
(2) 

Xm + X H + + X P H 2 + ZoH 2 + ^Hc + x c 

c s = \/fcBTj//imH, (3) 

where /ih and /1 are measured in units of the hydrogen mass 
7TiH- Note this this expression neglects the mass of elec- 
trons, and assumes that emitting species contribute negli- 
gibly to the mass. Two additional quantities are the gas 
specific heat at constant volume Cv,n and at constant pres- 
sure c Pi h, which for convenience we express per H nucleus 
rather than per unit mass or per unit volume. Thus c^h 
and c Pi h have units of energy over temperature, and can 
be converted to the usual values per unit mass simply by 
multiplying by a fiuma- Calculation of the specific heats 
requires some care when the chemical composition includes 
molecular hydrogen. I discuss this topic in detail, and derive 
DESPDTIC's expressions for c„,h and c Pi h, in Appendix 1X1 

The final quantity one can compute is the clumping 
factor / c i for the cloud, which represents an enhancement in 
the rates of all collisional processes due to non-uniformity 
of the gas. The quantity hh is the volume-averaged density 
over the cloud, but in a non-uniform cloud the density nn(x) 
at any position x may be higher or lower than this. Since the 
rate of collisions per unit volume at a given position varies 
as «h (a:) 2 , the rate of collisions per H atom in a non-uniform 
cloud exceeds that in a uniform cloud by a factor 



/cl = 



Wxf) 
n H 



(4) 



where the angle brackets indicate an average over the 
cloud volume; thus / c i is simply the factor by which the 
mass-weighted mean density exceeds the volume-weighed 
mean density. For a supe rsonically turbulen t medium, this 
factor is approximat e ly |Ostriker. Stone fc Gammiei 12001 



Padoan fc Nordlundl I2OO2I; also 



20081. iFederrath. Klessen fc Schmidt! 



Price. Federrath fc Bruntll201ll ) 



)01; 
one! 



Lemaster fc Ston 
<|2008l ). 



fcl&y/l + 0.75(T$, T /c2 



and 



(•>) 



2.2 Heating and Cooling Processes 

Given this physical model, one can proceed to calculate sev- 
eral routes by which the gas and grains can gain or lose en- 
ergy. In the following description, all heating, cooling, and 
energy exchange rates are given as energies per H nucleus 
per unit time. 



2.2.1 Gas 

2.2.1.1 Ionization Heating First, the gas can gain en- 
ergy through ionization heating; in this process primary elec- 
trons with energies produced when the gas is ionized by cos- 
mic rays or hard x-rays thermalize, adding energy. The rate 
at which this process adds energy is given by 



1 ion — £, Qion , 



(6) 



where q lon is the energy added per primary ionization. 
The value of qi on in turn depends on the bulk chemical 
composition of the gas, which determines how much of a 
primary electron's (a 37 eV of energy is lost via radiation 
rather than transformed into heat. T his problem has been 
discu ss ed by a number of author s ilDalgarno fc McCra' 



1972; Glassgold fc Langerl Il973l: IWolfire et al.l Il99 



ay 
15; 



Dalgarno. Yan fc Liul 11991 IWolfire.ll ollcnbac h fc McKeel 



201ol ; iGlassgold. Galli fc Padovanill2012l ). In predominantly 
atomic regions, the main pathway to thermalization is 
Coulomb scattering of the primary electron off other 
free electrons, and collisional excitation of H and He by 
the primary electron followed by collisional de-excitation 
of the excited atom. In this regime DESPOT IC uses the 
approximation recommended by iDrainel i|2011f) . 



gion.m « 6.5 eV + 26.4 eV f 

\X e + 



0.07 



1/2 



(7) 



In molecular regions the situation is far more compli- 
cated due to the additional thermalization channels provided 
by excitation of the rotational and vibrational levels of H2 
(followed by collisional de-excitation), by dissociation of H2, 
and by chemical heating, in which primary electrons produce 
reactive ions such as H^", H + , and He + that subsequently un- 
dergo exothermic reactions with neutrals such as CO, H2O, 
and O. In this case qi on becomes a complex function of the 
gas density and temperature, and the abundances of various 
species, and ranges from ~ 10 — 20 eV as these quantities 
change (jGlassgold. Galli fc Padovanill2012r ). Given the com- 
plexity of the problem, and the level of inaccuracy inherent 
in any one-zone model, DESPOT IC relies on a simple piecewise 
fit to the numerical results of IGlassgold. Galli fc Padovanil 
(|2012l . their Table 6) on the density-dependence of qi on in 
molecular regions: 



9ion,H 2 

eV 



10, 

10 + 3(log7i H -2)/2, 

13 + 4(logn H -4)/3, 

17+l(logn H -7)/3, 

18, 



log n H < 2 

2 < log n H < 4 
4^1og7i H <7 , (8) 
7 < log n H < 10 
logn-H ^ 10 



where the values of tih in the above expression are in units 
of cm -3 . 

To handle the case where the composition includes 
both molecular and atomic gas, DESPOTIC assumes that 
the atomic and molecular regions are physically separated 
(which, depending on the physical situation, may or may 
not be a good assumption). In this case the total heating 
rate can be computed simply by summing the heating rates 
in the atomic- and molecular-dominated regions, weighted 
by their number fractions: 



9ion = £HI?ion,HI + 2(x H 2 + £pH 2 )gion,H 2 



(9) 
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2.2.1.2 Photoelectric Heating Second, the gas can 
gain energy through grain photoelectric heating, whereby 
a primary electron ejected from a dust grain by a far- 
ultraviolet (FUV) photon thermalizes with the gas. Unlike 
cosmic rays, the FUV photons responsible for photoelectric 
heating can be attenuated by dust rather easily, and the 
photoelectric heating rate therefore depends on four factors: 
the strength of the ISRF, the abundance of dust grains, 
the amount of dust shielding, and the energy yield per 
photoelectron; as with cosmic ray heating, t he latter value 
has been estimated by numerous authors ilWat son 1972; 
de Jond[l977l; iTielens fc Hollenbacblll985l ; [Bakes fc Tielensl 



1994 



Wolfire et al.1 



2003) 



To account for dust shielding, 
which obviously varies from point to point within a real 
cloud, DESPOTIC uses the simple ap proximation proposed by 
iKrumholz. Lerov fc McKeei (|201ll ), whereby the 8 - 13.6 eV 
photons responsible for photoelectron production are con- 
sidered to be attenuated by half the mean extinction of the 
cloud. Since the dust opacity is relatively flat across this en - 
ergy range (~ 50% variation in the models of|Draine]|2003), 
we can assign a single cross section ct^pe, which is ~ 10 -21 
cm 2 H _1 for Milky Way dust. Thi s value is near the middle 
of the range found in the models of lDrainei (|2003l ). With this 
approximation, the photoelectric heating rate becomes 



r PE = 4.0 x 10" 26 x^e 



-(l/2)iV H <T e 



erg s 



H/ 



(10) 



2.2.1.3 Gravitational Heating A third possible source 
of heating is adiabatic compression. This obviously depends 
on the hydrodynamics of the flow, something that is not 
naturally included in a one-zone model like that used in 
DESPOTIC. However, this effect is calculable in the special 
case of compression due to gravitational contraction, as in 
protostellar cores for example. In this case the heating rate 
may be computed using the approxim ation introduced by 



may be computed using the approxim a 
iMasunaea. Mivama fc Inutsukal IJ1998T ). 



Ci c 2 /x H m H i/47rGp, 



(11) 



where C\ is a dimensionless constant of order unity that de- 
pends on the nature of the gravitational collapse. From their 
numer ical calculations, iMasunaga. Mivama fc Inutsukal 
(1998) find C\ ~ 1.0. Since in general most interstellar 
clouds are not in a state of collapse, by default DESPOTIC 
does not include gravitational contraction heating, and sets 
C\ = 0. However, users do have the option of overriding 
this default. 

2.2.1.4 Line Cooling The primary cooling mechanism 
for gas is line radiation. For each emitter species s, there is 
a rate of line cooling A s , so that the total line cooling rate 

is 



Alin 



E 



As 



(12) 



I defer a calculation of A s to § 12.31 



2.2.1.5 Dust- Gas Energy Exchange Finally, gas can 
either heat or cool by exchanging energy with the dust via 
collisions. The gas-dust energy exchange rate is given by 



* g d = tt g dfcinnTg' 2 (T d — T g 



(13) 



where a g( j is the grain-gas coupling coefficient and the sign 
convention is that positive values correspond to heating of 



the gas and cooling of the dust. Note the presence of the 
clumping factor / c i, since this is a collisional process. The 
coupling constant depends on the grain abundance, chem- 
ical composition, size distrib u tion, and charge state. For 
Milky Way dust, iGoldsmithl (J2001T ) recommends a value 
a g d = 3.2 x 10~ 34 erg cm 3 K~ 3 ' 2 for H 2-dominated regions, 
and IKrumholz. Lerov fc McKeei ( J2011T ) estimate a value of 
1.0 x 10~ 33 erg cm 3 K _3//2 for H i-dominated ones, with the 
difference arising due to the change in both the number and 
mean mass of free particles between H I and H2-dominated 
regions. 

2.2.1.6 Summary of Gas Heating, Cooling, and 
Thermal Evolution In summary, the time rate of change 
of the gas energy per H nucleus e 9 , sp is given by 

de 

-Jp = Tion + TPE + r grav - Al inc + *g d . (14) 

The corresponding time rate of change of the temperature 

is 



dT^ 
dt 



(Cv,H, C p ,h) 



(ri 0n + Tpe + r g 



A linc + * gd ), (15) 



where c„,h is the gas specific heat per H nucleus at con- 
stant volume and c Pi h is the specific heat per H nucleus at 
constant pressure. The parentheses indicate that one can 
use either Cu,H or c Pi h in the above equation, depending on 
whether one wishes to consider gas cooling isochorically or 
isobarically. 

2.2.2 Dust 

2.2.2.1 Cooling by Thermal Radiation Dust grains 
can lose energy via thermal continuum radiation. To com- 
pute the cooling rate, consider a population of spherical 
grains with distribution of radii a g given by dn/da g , where 
we normalize the distribution function such that rid — 
J(dn/da g ) da g is the total number density of dust grains. 
Let Q v (a g ) be the absorption efficiency for absorption of 
radiation of frequency v, so that the cross section of the 
grain to radiation of frequency v is a„(a g ) = Tra^Q(u). Fur- 
ther let (Q(a g )) T = J B v (T)Q v (a g ) dv / ' / B V (T) dv be the 
Planck- weighted mean efficiency, where B V (T) is the Planck 
function evaluated at temperature T. Given this definition, 
we can write the rate of thermal radiation cooling from dust 
grains of temperature Td as 



A rf 



1 f dn , . .. 2 , 

/ ~, — {Qv(a g ))TTra da g 

nH J da g 

a d (T d )caT d , 



caTd 



(16) 
(17) 



we have defined the term in square brackets to be the mean 
dust cross section per H nucleus (Jd(Td). This expression as- 
sumes that the cloud is optically thin to its own cooling 
radiation; we treat the optically thick regime below. We ap- 
proximate that Od{Td) will vary as a powerlaw with Td, and 
we therefore write 

a d (T d ) = cr d ,io (j^) . (18) 

For Milky Way dust, typical opacities are ad 10 ~ 2 x 10 _ 
cm 2 H _1 l|Pollack et al.lll994lSemenov et al.ll2003l ). and for 
temperatures T d such at hc/{k B T d ) = 0.14(T d /10 K)" 1 cm 
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is much larger than the typical grain size, we expect fi — 2; 
detailed grain models show that thi s expectation holds up 
to T d « 150 K (jSemenov et al.ll2003h . DESPOTIC leaves both 
(Td,io and P as user-settable parameters. A naive expectation 
is that, at sub-Solar metallicities, o"d,io oc Z' d , where Z' d is 
the dust abundance relative to Solar. 

The above estimate is valid only as long as the cloud is 
optically thin to its own cooling radiation, which is true only 
as long as Od,io(Td/10 K) Ah < 1. Given the small value of 
o"d,io for Milky Way dust, departures from the optically thin 
regime do not begin until extremely high column densities. 
However, there are circumstances, for example in the molec- 
ular clouds of starburst galaxies, where T d and Ah can be 
high enough to render the optical depth to cooling radiation 
large. A truly accurate calculation of the cooling rate in 
this regime requires a multi-zone numerica l treatment with 
a rad iative transfer code such as dusty l|lvezic fc Elitzurl 
1 19971 ) or SteinRay (jSteinacker et alj 2003). or a sophisti- 
cated analytic approximation (e.g. IChakrabarti fc McKed 
120051 ). However, we can obtain a very crude treatment of the 
optically thick regime by noting that the maximum possible 
cooling rate for the cloud is simply irR 2 caT d , the blackbody 
rate for a sphere of radius R = (3/4)A r H/«-H equal to the 
cloud radius. Rewriting this as a rate per H nucleus, the 
maximum possible dust cooling rate is 



tion, the heating rate of grains due to the ISRF is 



A^ thick — 



coTl 



DESPOTIC adopts the approximation 
A d = min(A dithin , A d , thick ). 



(19) 



(20) 



2.2.2.2 ISRF Heating Grains can be heated by ab- 
sorbing the interstellar radiation field produced by stars. 
To compute the rate of dust heating from the ISRF, 
we must perform a calculation similar to that for A d . 
In analogy to (Q„(a g )) T , we define (Q„(a 9 ))rsRF = 
j u v ^ishf Q v {a g ) du / j u v ^ishf du, where u„,isrf is the en- 
ergy density of the ISRF at frequency v, as the ISRF- 
averaged absorption efficiency. In general (Q„(a 9 ))isRF 3> 
(Qi/(« s ))t. Thus, unlike in the case of thermal cooling where 
optical depth effects are important only in extreme circum- 
stances, attenuation of the ISRF will be important even 
at modest column densities. As with photoelectric heating, 
it is clear that there is no single value that describes the 
rate of dust heating within an optically thick cloud; heating 
rates will be high at the edge and low at the center. More- 
over, unlike in the case of photoelectric heating, the range 
of photon energies responsible for heating is quite broad, 
with half the heating coming from photons with wavelengths 
> 0.31 /im even for the unattenuated ISRF (B. Draine, 2013, 
priv. comm.). As a result, the spectrum of the heating field 
changes as one moves into a cloud and shorter wavelength 
photons are selectively attenuated. Consequently, in addi- 
tion to the geometric uncertainty, there is an additional one 
in the choice of dust cross section to assign. In order to 
maintain simplicity, DESPOTIC does not attempt to treat this 
problem in detail, but instead uses the same approximation 
as for photoelectric heating, i.e. that the characteristic heat- 
ing rate is to be computed assuming an attenuation equal to 
half the mean value for the cloud, using a single grain cross 
section to compute the attenuation. With this approxima- 



Tisf 



1 
"II 



dn 
da„ 



{Qiy(a g )}isnFTva g da s 



' c-uisRFe 

24 



-^dJSRF-^H/ 2 



= 3.9 x W~ M xZ' d e 



~' T d,ISRF N H/ 2 



(21) 
ergs" 1 H~ 1 ,(22) 



where Z' d is the dust abundance relative to the Milky Way 
value, misrf = X u mw is the energy density of the ISRF, 
umw is the energy density for the Milky Way's ISRF, ctisrf 
is the cross section we assign for IS RF attenuat i on, an d the 
numerical coefficient is taken from iGoldsmithl ([2001). The 
choice of Uisrf is somewhat difficult for the reasons stated 
above, and if very high accuracy is desired it should be com- 
puted on a case-by-case basis. However, a reasonable default 
for Milky Way dust is ctisrf = 3 x 10 -22 cm 2 H~\ which 
is roughly halfway between the values appropriate for the 
unextincted ISRF and the value expected for an ISRF ex- 
tincted by an optical depth of 2 in V band (B. Draine, 2013, 
priv. comm.). 

It is worth noting that, because the ISRF is exponen- 
tially attenuated by dust, when <t £ j,isrfA'h > lwe are likely 
to find that Fisrf is negligibly small even when \ i s very 
large. In this circumstance, the ISRF is so thoroughly atten- 
uated that none of it reaches the cloud interior where we are 
computing the temperature. However, if this happens, the 
hot outer parts of the cloud that are directly exposed to the 
ISRF will heat up and generate a background infrared field 
within the cloud interior. If the cloud is optically thin to IR 
cooling radiation the intensity of this field will be low and 
it can be neglected as a heat source. If the cloud is optically 
thick to IR, on the other hand, the background IR field will 
build up, and will heat the cloud interior. DESPOTIC pro- 
vides a mechanism to handle this phenomenon by including 
an infrared radiation field (see the following section), and in 
circumstances where ISRF heating is negligible, heating by 
the infrared radiation field should take its place. As for the 
case of the cooling rate when the cloud is optically thick to 
IR, calculating the intensity of the background field in this 
circumstance requires a more sophisticated model than the 
one-zone treatment that DESPOTIC provides. However, we 
can solve the limiting case of an extremely optically thick 
cloud subject to external heating. If such a cloud absorbs 
all of the background ISRF incident on its surface, the to- 
tal heating rate is 7t.R 2 ciiisrf, and the heating rate per H 
nucleus is 



Tisf 



CUISRF 

Ah 
5.3 x 10" 



Nn 



erg s H 



(23) 



Equating this with the limiting cooling rate for an extremely 
opaque cloud, A^ thick, gives an equilibrium temperature for 
both the dust and the infrared radiation field 



Td.thic 



T IB 



(*) 



1/4 



2.1 X 1/4 K 



(24) 



i.e. the dust and IR radiation field within the cloud reach a 
temperature such that the radiation energy density within 
the cloud is equal to the ISRF energy density outside it, as 
expected for a blackbody. 
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2.2.2.3 Heating by Infrared Radiation and the 

CMB The final source of radiative energy for dust is the 
background thermal radiation field, and the CMB. Since 
both of these sources of radiation are thermal, they may 
be handled using exactly the same mechanics as thermal 
radiative cooling. The heating rate is therefore 



10 K 



fT Tl 

TdJR = G"d, 10 

ry _ ( TCMB \ 

1 d.CMB — O"d,10 I ..» t7 j 



P 



caT r . 



P 



caT% 



(25) 
(26) 



2.2.2.4 Line Heating In addition to emission and ab- 
sorption of continuum radiation, there are two additional 
processes that can heat and cool dust grains. The first 
of these, collisional exchange with the gas, is discussed in 
§ 12.2.11 The other is absorption of line photons emitted by 
the gas. If we let 



dn 2 , 
Cd,u = — / -, — Ka„ 

7lH / «<2e 



, da a 



(27) 



be the population-averaged grain cross section per H nucleus 
at frequency v, then the mean optical depth of the cloud 
to line photons at frequency v is Td,u = Nno~d,v In prin- 
ciple one could use a detailed grain model to obtain aa,v 
at the frequencies of all the relevant lines. However, this 
procedure would be cumbersome, and is likely unimportant 
for most clouds since, not surprisingly, both cooling radia- 
tion and observable emission tend to be dominated by lines 
at frequencies such at clouds are optically thin. Nonethe- 
less, to approximate the effects of clouds becoming opti- 
cally thick to line radiation, DESPOTIC approximates aa,v 
by °d,io (v/209, GHz) , where v is the line frequency, and 
208 GHz is (k B /h) multiplied by 10 K. With this approxi- 
mation, and using the same expression for the line photon 
escape probability versus optical depth as discussed below 
in § 12.31 we obtain the final heating rate of the dust due to 
absorption of line photons: 



-L d, line — / V Pd.s.ij ) *ys,ij 

Pd,s,ij 



1 + fiVH0- d ,io Ky/208 GHzr 



(28) 
(29) 



where v Si ij is the frequency of the line produced by atoms 
/ molecules of species s transitioning between states i and 
j (see § I2.3[l . Pd,s,ij is the escape probability for a photon 
corresponding to line ij computed using the dust optical 
depth, and the sum runs over all species s and level pairs ij. 



2.2.2.5 Summary of Dust Heating and Cooling 
Processes Collecting all the terms for dust, the total rate 
of change of the dust specific energy per H nucleus is 



d(:, 



d,sp 



(It 



= Tisrf + r<j,] 



IVcmb + r<j,iR — Ad — * g d- (30) 



In principle one could consider time-dependent temperature 
evolution of the dust as well as of the gas, but since the 
specific heat of the dust is far less than that of the gas, 
and is a complex function of the properties of the grains, 
DESPOTIC does not treat this case. Instead, it assumes that 
the grain population is always in thermal equilibrium. 



2.3 Level Populations and Line Radiation 

2.3.1 Level Populations in Optically Thin Clouds 

Calculating the line cooling rate requires determining the 
level populations for all emitting species. Consider an emit- 
ting species s, and let Ei be the energy of the ith quantum 
state of that species, where the states are numbered by en- 
ergy so that Ei < Ei+i for all states i. The degeneracy of 
state i is gi, and the Einstein coefficient describing the rate 
of spontaneous radiative transitions from state i to state j 
is Aij, where Aij = for i ^ j. Finally, let k Py ij be the rate 
coefficient for collisional transitions from state i to state j 
induced by collisions with some collision partner p; the up- 
ward and downward rate coefficients obey the usual rela- 
tionship kpji = {gi/gj)k p ,ij exp(—AEij/kT g ), where i > j 
and AEij = \Ei — Ej\. By convention k Pt ij — for i = j. 

For our species of interest, we wish to solve for the frac- 
tion fi of atoms / molecules in state i, when that species is 
mixed with a gas of a given bulk composition, number den- 
sity riH, and gas temperature T g , and the cloud is immersed 
in a sea of cosmic microwave background photons. If the 
cloud is optically thin to photons at the frequencies of the 
lines connecting the various states, in statistical equilibrium 
the various level populations are determined implicitly by 
the conditions that the rate of transitions into and out of 
each level balance: 



£/< 



Qji "T" I J- "T" '"'Y, 3^)-"- ji ' <tv»~f ,ij ■**- i ij 

9i 



*£ 



where 



Qik + (1 + n 1: ik)Aik H n liki A ki 

9i 



1 



7,lJ ' exp(A£ ij /fc s r C MB) - 1 

Qij — Jcl ^H / £p fcp , ij 



(31) 

(32) 
(33) 



are the photon occupation number at the frequency of the 
line connecting states i and jp and the rate of collisional 



3 Naively one would think that, in a cloud that builds up a signif- 
icant trapped infrared radiation field, then the photon occupation 
number should also include a contribution from this field, of the 
same form as equation l |32jl but with Tcmb replaced by T ra( i dust- 
However, this is often not the case, for the following reason. Even 
in high column density environments where a significant dust- 
trapped infrared radiation field builds up, the spectrum of this 
radiation field is often not Planckian at low frequencies. This is 
because the dust opacity generally falls as u 2 at low frequencies, 
and so even if the dust is opaque to radiation near the peak of 
the spectral energy distribution, it is usually transparent at low 
frequencies. This results in a radiation spectrum that is Planck- 
ian at higher frequencies but very sub-Planckian at low frequen- 
cies, and thus has a much lower photon occupation number that 
a true blackbody like the CMB. A fully accurate calculation of 
level populations would account for this effect by solving for the 
frequency-dependent dust-mediated radiation field and using the 
appropriate photon occupation number to calculate the level pop- 
ulations. However, as noted above, it is not feasible to determine 
the dust radiation field accurately in a one-zone model. I there- 
fore choose to optimize the accuracy of DESPOTIC for the case of 
lines at frequencies where the dust is optically thin, since these 
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transitions between states i and j summed over all collision 
partners p. Here x p is the abundance of a given collision 
partner relative to tih, and the collision partners considered 
by DESPOTIC are H, He, pH 2 , oH 2 , e, and H + . As usual, col- 
lision rates are multiplied by the clumping factor / c i. The 
left-hand sides of equations (|3ip describe the rate of transi- 
tions into state i from all other states j, with the first term 
representing the rate of collisional transitions, the second 
representing the rate of radiative transitions (including both 
spontaneous and stimulated emission), and the third term 
describing the rate of absorptions. The right-hand sides rep- 
resent the rate of transitions from state i to all other states 
k, with the three terms again representing collisional transi- 
tions, spontaneous and stimulated emission, and absorption. 
These equations are supplemented by the constraint equa- 
tion 



E* 



i, 



(34) 



and together equations (|3ip and (|34p constitute a complete 
system. 

2.3.2 Level Populations in Optically Thick Clouds 

If the cloud is optically thick, these equations must be mod- 
ified to account for the fact the effects of the trapped ra- 
diation field that builds up inside the cloud. To handle this 
case, DESPOTIC uses the standard escape probability approx- 
imation, in which the level populations are assumed to be 
uniform, and every transition ij is assigned an escape proba- 
bility /3ij, which gives the volume-averaged probability that, 
when an atom or molecule radiatively decays from state i 
to state j, the associated photon will escape from the cloud 
rather than being resonantly absorbed within it. With this 
approximat i on , the modified equations simply become (e.g. 
iDraindlJOlTI ) 






Qii \ Pji\*- ~r Tl-y ji )/±ji ~r Pij Tl^i ij/\%j 

9j 



q%k + /3ife (1 + n~/,ik)Aik + fiki H 7l fcjAfcj 

9i 



(35) 



The escape probability may be computed using several 
possible approximations, which are appropriate for different 
cloud geometrie s. By d e fault, DESPOTIC uses the approxi- 
mate result from lDraind (|201ll ) for uniform spherical clouds, 



Pa 



<)< 



a'v 

AijXij 



gj 4(27r) 3 / 2 atot 



XsNnfj 1 



fidj 
fj9i 



(36) 



(37) 



where Tij is the optical depth corresponding to a column 
A^h, Aij = hc/Aij is the wavelength of transition ij, a to t — 
\/o- 2 NT + Cs/pLs, /^s is the molecular weight of the emitting 



are, obviously, the lines that are most important for both cool- 
ing and observation. This choice dictates that the dust radiation 
field be ignored when computing the level populations, on the 
basis that its photon occupation number will be small. However, 
this choice does limit the accuracy of DESPOTIC for lines where 
infrared pumping is important, as discussed in more detail in §[5] 



species in units of thh, and x 3 is the abundance of the emit- 
ting species per H nucleus. Note that, in the expression for 
Pij, the coe fficient on rg differs by a f actor of (3/4) from 
that given in iDrainel i|201 if) because iDrainei defines Tij using 
the center-to-edge rather than the mean column density. 

The code can also use one of two other approxi- 
mations. For a slab geometry, the escape probability is 
l|de Jong. Dalgarno fc ChuHl975l l 



A 



-3-ry 



ij ,slab 



37V, 



(38) 



but now iVn is 
interpreted as the column density of the slab rather than 
the mean column density of a sphere. Finally, DESPOTIC 
can use the large velocity gradient (LVG) approximation, in 
which the escape probability is computed from a Sobolev 
approximation and the geometry is therefore irrelevant 
(jde Jong. Boland fc Dalgarndl 19801 ). In this case 



Pij,ISVG 



Tij ,LVG 



1 _ g _T ij,LVG 

Tij, LVG 

9i A a x g M , (. hg 3 
gj 8tt \dv r /dr\ V fcgt 



(39) 

(40) 



For whichever choice of geometry, the above equations 
determine Pij in terms of fi and other known quantities, 
and together with equations p4[) and (|35|l they again form 
a complete system that may be solved for fi. In the optically 
thin limit, /3ij — > 1 for all ij, and equations ([35jl reduce to 
equations (|31|l . 



2.3.3 Line Cooling Rates 

Given a set of level populations determined by solving the 
equations given in the previous section, the cooling rate of 
the cloud due to emission by line ij of species s is given by 



l*-s,ij — Pi' 



(1 + n,~ f ,ij)fi "n>i,ij)j 

9j 



ij ^-E^ijX s 



(41) 



Note that this is the net cooling rate, in that the first term 
in brackets represents the rate of spontaneous plus stimu- 
lated emission per emitting atom / molecule, while the sec- 
ond term is the rate of absorption of background photons. 
Thus A S: ij is the rate of energy loss via line emission minus 
the rate of energy gain from absorption of the background 
radiation field. If T g < Tcmb, then A 3 ,ij will be negative, 
indicating a net gain in energy. The total cooling rate A 3 for 
species s is simply the sum over all level pairs, 



As 



E A * 



(42) 



2.4 Line Shapes 

DESPOTIC's final major capability is calculating the profiles 
of spectral lines. In general this is not a useful calculation in 
a one-zone escape probability model; since the level popula- 
tions in such a model are assumed to be uniform, the result 
is necessarily rather uninteresting, and is simply given by 
the usual solution to the transfer equation for media with 
enrission and absorption coefficients that are independent of 
position. However, one can relax the assumption of uniform 
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-VR 2 - d 2 




s = Vi? 2 - d? 

To observer 



Figure 1. Diagram of the geometry used by DESPOTIC when 
calculating line shapes. The circle shows the cloud, with radius 
R = 3-/Vu/4nH, and the dashed line is the observer's line of sight 
through it. 



level populations by making another one: that the species is 
in LTE, and that the temperature T is a known function of 
position|j Solving for the shapes of lines in this limit allows 
the code to compute pCygni and inverse pCygni profiles, 
among other applications. 

Consider a spherical cloud following DESPDTIC's general 
model, and consider a line of sight passing through it at 
an offset distance d from the cloud center (see Figure [T]). 
We let n a , T, v, and ctnt be the number density of the 
emitting species s, gas temperature, the radial velocity, and 
the non-thermal velocity dispersion, and in general all can 
be functions of r. Now consider a spectral line of this species 
connecting an upper state u to a lower state £. Under the 
assumption of LTE, number densities of species s in the 
upper and lower states are 



" t: 



gee 



-E t /k B T 



Que 



-E u /k B T 



Zs(T) 



Zs{T) 



(43) 



where gt is the degeneracy of state i, Ei is the energy of the 
state, and Z a (T) is the partition function for species s at 
temperature T. 

The equation of radiative transfer along the chosen line 
of sight reads 



dh 
ds 



> 



K V -t V 7 



(44) 



where s is the position along the line of sight, defined such 
that s = is the cloud midplane, and integration through 
the cloud proceeds from s — — \/R 2 — d? to +\ / 'R 2 — d? (Fig- 
ure!]]). The emission and absorption coefficients are given by 



9u A 

= — nt—A ul (, 

gi 07T 

= k v B v (T) 



(45) 
(46) 



In principle one in fact needs to know only the excitation tem- 
perature T ex for the two levels that produce the line, together 
with the number density the atoms or molecules that are in the 
lower state rip. However, in practice it is unlikely that one will si- 
multaneously know T ex and np in any situation other than when 
the levels are in LTE, and thus I limit the discussion to this case. 
If one does in fact know T ex and ni, it is trivial to perform a 
calculation for that case simply by setting the level populations 
to their LTE values at T ex , and adjusting the overall density of 
the species so that ni has the desired value. 



where A = hc/AE is the wavelength of the transition, AE — 
E u — Ei is the energy difference between the levels, A u £ is 
the Einstein coefficient for the transition, and </> v is the line 
profile. This is given by 



V2t< 
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(47) 



where o v and vq are the dispersion in frequency and the 
frequency of line center, given by 
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(48) 
(49) 



where /j, a is the mass of a particle of species s, measured in H 
masses. The transfer equation may be non-dimensionalized 
via the change of variables. We let x — s/R be the dimen- 
sionless position, / = i//(AE/h) be the dimensionless fre- 
quency, Xf — I v /(A u en a (R)hR) be the dimensionless in- 
tensity, and we normalize all the position-dependent quan- 
tities to their values at the cloud edge: n' s = n s /n s (R), 
t = T/T(R), ip - o-nt/ctnt(-R), and u = v/v(R). With these 
definitions, after some manipulation the transfer equation 
becomes 
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and we have defined the dimensionless ratios 
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(50) 

(51) 
(52) 
(53) 

(54) 
(55) 
(56) 



The dimensionless intensity 1/ at any dimensionless fre- 
quency / may be obtained by integrating equation ()50[) from 
x = -y/l - (d/R) 2 to x = +\/l - (d/R) 2 subject to the 
boundary condition Xf — B v (Tcmb) / (A u tn s (R)hR) at the 
lower limit of integration. 



3 CODE ARCHITECTURE AND 
ALGORITHMS 

In this section I describe the architecture of the DESPOTIC 
code and the algorithms it uses to solve the equations intro- 
duced in the previous section. 

3.1 Overall Architecture 

DESPOTIC is a library intended not only to be used for stand- 
alone calculations, but also to allow easy extensibility, easy 
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integration with other codes, and to allow users to conduct 
interactive, exploratory calculations. To this end, DESPOTIC 
is implemented as a Python package. This package defines 
a series of classes, which store physical and chemical prop- 
erties of clouds, and which define a series of methods for 
calculating quantities of interest from those properties. A 
full listing of the code's methods and procedures is given in 
the User Guide that is included with the DESPOTIC package. 
This architecture also permits a very high level of abstrac- 
tion, such that many useful computations can be performed 
with no more than a single line of code on the part of the 
user. I provide examples of such computations in § [3] Im- 
plementation in Python also provides access to a number 
of other useful features; for example, the cloud descriptor 
classes that DESPOTIC defines can be pickled via the usual 
Python tools, making it possible to save intermediate states 
of computations in a portable, machine-independent format. 
To achieve high performance, DESPOTIC makes 
extensive use of the ability of the numPy and sciPy 
libraries to interface with the fast, optimiz ed numeri- 
cal libraries LAPAClfl ilAnderson et all 1 19991 ). MINPACf 
More. Garbow fc Hillstroml Il980l ) 



and 



ODEPACf 



Hindmarshl ll983T ). All the computationally- intensive 



operations performed by DESPOTIC are handled through 
such interfaces. While an interactive Python code like 
DESPOTIC will never be as fast as a code written in a lower- 
level language like Fortran and compiled with an optimizing 
compiler, the use of these interfaces erases enough of the 
difference to make DESPOTIC computations reasonably 
rapid. It is hard to provide a quantitative estimate of code 
execution times for DESPOTIC routines, since as I discuss 
below the most computationally-intensive ones require 
iterative methods, and the time required for such a solution 
is a strong function of the quality of the starting guess. 
Nonetheless, I give a general idea of code execution times, 
as tested on a single processor of a modern workstation, 
for some example applications in § [4] Individual instances 
of DESPOTIC classes use internal private storage, and thus 
are thread-safe should a user desire to use threading to 
accelerate the calculation of large grids of models via the 
standard Python threading interface. Threading of internal 
DESPOTIC calculations for single clouds will be added in a 
future release. 



3.2 Capabilities and Algorithms 

3.2.1 Optically Thin Level Populations and Line 
Luminosities 

The most basic capability of DESPOTIC is to compute level 
populations and line luminosities for an emitting species 
embedded in a cloud of specified physical properties (tih, 
T g , ctnt, abundances, etc.) - thi s is the same computa - 
tion performe d by codes li ke RADEX (van der Tak et alj|2007l ) 
and lineLum (Krumholz & Thompson 2007), and the lat- 
ter is the direct ancestor of the corresponding portion of 
DESPOTIC. Appendix [B] provides a direct comparison be- 
tween DESPOTIC and RADEX. In the optically thin regime, 



this operation requires that the code solve the system of 
equations ([31) and ((34} . For computational purposes it is 
convenient to rewrite the system as a matrix equation. Con- 
sider a species s for which we track TV distinct energy levels. 
With some manipulation, equations (|31|l and (|34[) may be 
rewritten aqfj 

Mf = b (57) 

where M is an (TV + 1) X TV matrix whose elements are 

My = — 5ij + <5i,jv+i 



+ 



J2k [l ik + (^ + n -t,ik)Aik + ^n 7i fciA fci ] 
b is a vector of length TV + 1 whose elements are 



*i. 



JV+l 



(58) 



(59) 



and / is a vector of length TV whose elements are the frac- 
tional level populations fi. By convention g^ = Aij — 
n i,ij — f° r * = N + 1 or j — TV + 1. The off-diagonal 
elements of matrix M in rows i ^ TV have a simple physi- 
cal meaning: element ij is the rate coefficient for transitions 
(adding both radiative and collisional processes) into state i 
from state j, normalized by the sum of the rate coefficients 
for all transitions out of state i to any other state. The final 
row of M, i — TV + 1, implements the constraint equation 
that the sum of all fractional level populations is unity. 

DESPOTIC first constructs the matrix M from the spec- 
ified cloud properties. In most cases, it then immediately 
solves equation (|57|l using the LAPACK routine lstsq. How- 
ever, robust numerical solution of equation (|57|l requires con- 
siderable care, because when the transition probabilities into 
certain states are very low, the matrix M can by extremely 
ill-conditioned, making accurate numerical solution impos- 
sible. To avoid this problem, DESPOTIC evaluates the condi- 
tion number of M, and, if it is excessively large, it employs a 
number of techniques to reduce it before proceeding to nu- 
merical solution. Those techniques are detailed in Appendix 

El 

Once the level populations have been found, it is triv- 
ial to compute the radiation rate per H nucleus from each 
line using equation (|4ip with f3ij = 1. In addition to the 
line luminosity per H, DESPOTIC also computes the emergent 
frequency-integrated intensity I 3t ij and velocity-integrated 



ft 



/^A^TVh 

47T 



Tb,. 



h/kl 



' ln[l + 2/u^/c 2 / s 



(60) 

(61) 



5 http://www.netlib.org/lapack/ 

6 http://www.netlib.org/minpack/ 

7 https : //computation . llnl . gov/casc/odepack/odepack_home . html 



Note the factor fid,s,ij in equation (|60[) . which accounts for 
absorption of line radiation by dust internal to the cloud. 
This effect need not be included when calculating the level 



8 Note that DESPOTIC does not use the standard procedure in 
the stellar atmospheres community of recasting the equations in 
terms of departure coefficients. This choice is motivated by the fact 
that, for most of the calculations for which DESPOTIC is intended, 
most of the states of most species will be very far from LTE. This 
vitiates any advantage to recasting the equations in departure 
coefficient form. 



© 0000 RAS, MNRAS 000, 000-000 



10 Krumholz 



populations, under the assumption that any line photon ab- 
sorbed by dust will not be re-emitted in resonance with the 
line and thus cannot cause an absorption elsewhere in the 
cloud. This assumption is well-justified for the infrared and 
radio lines for which DESPOTIC is specialized, since absorp- 
tion opacities exceed scattering opacities at these frequen- 
cies by many orders of magnitude. However, dust absorption 
must still be included when calculating the observable in- 
tensity emerging from the cloud, since photons absorbed by 
dust will ultimately be emitted as thermal continuum rather 
than lines. 



3.2.2 Optically Thick Level Populations and Line 
Luminosities 

DESPOTIC can also calculate level populations and line lumi- 
nosities for optically thick clouds, which requires solving the 
system of equations (|34[) and (|35[) . together with equation 
(|36p . (|38p . or (|39p to approximate the escape probability. For 
specified escape probabilities /3y , the procedure for calculat- 
ing the level populations is identical to that given in § 13.2.11 
for optically thin clouds, except that matrix M becomes 

Mij — —5ij + 5i t N+i 



+ 



(62) 



However, the escape probabilities fiij are not known in ad- 
vance, and they and the level populations must instead be 
computed iteratively. Again, by convention, /% = for 
i = N+l. 

Let fl be the current best guess for the level pop- 
ulations after n iterations. At every step of the iteration, 
DESPOTIC uses /} to compute a new estimate for the es- 
cape probabilities p£ of all lines from equation (|36p . (138|) . 
or (|39p . constructs the matrix M following equation (|62[l us- 
ing pQ' , and then solves equation (f57|) to obtain a new set of 
level population estimate fl ■'. DESPOTIC then checks if the 
level populations have converged by computing the absolute 
and relative residuals 
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(63) 
(64) 



and comparing them to specified tolerances. If the residuals 
exceed the specified tolerances, DESPOTIC generates a new 
set of level populations 



fl n+1) =Df^ + (l-D)ri n \ 



(65) 



where the damping factor D is in the range (0,1]. Larger 
values of D represent more aggressive attempts to con- 
verge to the solution rapidly, at the cost of a higher risk 
of non-convergence. DESPOTIC chooses a default D = 0.5, 
but this value can be altered by the user, and in calculations 
where level populations are likely to be computed repeatedly 
(for example when computing thermal equilibria), DESPOTIC 
catches non-convergences automatically and attempts to re- 
compute using a smaller value of D, thereby preventing the 
entire computation from being derailed. 



to either the currently-stored level populations for a given 
cloud or, if none are available, their LTE values for the gas 
temperature. Initializing to the currently-stored level popu- 
lations ensures that, when level populations must be com- 
puted repeatedly under physical conditions that very only 
slightly, as in many of the examples given in § 3] the ini- 
tial guess will be close to the correct level populations, and 
convergence will be rapid. 

Once a converged solution for the level populations is 
found, DESPOTIC calculates the line luminosities using equa- 
tion (|4ip . and the integrated intensity and brightness tem- 
perature emerging from the cloud via equations (|60p and 
(|6ip . This is identical to the optically thin case, except that 
the escape probabilities f3ij may not be unity. 



3.2.3 Cooling Rates, Thermal Equilibria, and 
Time- Dependent Temperature Evolution 

In addition to computing line luminosities and level popu- 
lations, DESPOTIC can also compute the heating and cooling 
rates of gas and dust. It does so by evaluating all the terms 
in equations (|14p and (|30| l: since one of these terms is Aii nc , 
this procedure entails solving for the level populations and 
escape probabilities as described in § 13.2.21 

DESPOTIC can also solve for equilibrium dust and gas 
temperatures. DESPOTIC obtains these values by setting 
de g , sp /dt — in equation (|14p and ded, sp /dt = in equa- 



To start the process, DESPOTIC initializes by setting f 



(0) 



tion (|30p . The user can also add arbitrary additional heating 
and cooling terms to either equation, to represent processes 
not modeled by DESPOTIC (e.g. endothermic or exothermic 
chemical reactions). At the discretion of the user, DESPOTIC 
can fix either T g or Td and solve for the other, or it can 
solve for both simultaneously. If either T g or Td is fixed, 
DESPOTIC solves the equations using the secant method. If 
neither is fixed, it solves for T g and T<j simultaneously using 
the MINPACK routine hybrdl, which implements the Powell 
hybrid method. 

Finally, DESPOTIC can compute the time-dependent 
thermal evolution of a cloud. Starting from an initial gas 
and dust temperature, DESPOTIC can integrate equation (|15[1 
for the gas temperature evolution. At the user's discretion, 
the calculation can be done either isochorically or isobari- 
cally. When evaluating the heating and cooling terms that 
appear on the right-hand side of equation (|15p . DESPOTIC 
assumes that both the level populations and the dust tem- 
perature reaches equilibrium instantaneously; the former are 
computed via the procedure described in § 13.2.21 and the 
latter by the solution to equation (|30[) with ded, sp /dt = 0. 
DESPOTIC also calculates the temperature-dependent specific 
heat c^h or c Pi h on the right-hand side using equation (IA9[) . 
It integrates equation (|15p using the 0DEPACK routine lsoda, 
which automatically evaluates the stiffness of the system, 
and solves using a predictor-corrector method for non-stiff 
problems and backward differentiation formula methods for 
stiff problems. 



3.2.4 Line Profiles 

DESPOTIC's final major capability is calculating line profiles 
for species in LTE. When performing this calculation, it ac- 
cepts user-specified profiles for the number density of the 
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emitting species, the bulk velocity, the non-thermal veloc- 
ity dispersion, and the temperature as a function of radius. 
From these inputs, plus the identity of the line whose pro- 
file is to be computed, it calculates all the dimensionless 
quantities given in equations ((54)) - (f56|) . and then numer- 
ically integrates equation (|50|) at a range of user-specified 
frequencies or velocities. The integration is performed via a 
call to the ODEPACK routine lsoda. DESPOTIC then returns 
the CMB-subtracted intensity and brightness temperature 
as a function of frequency / velocity. 



3.3 Atomic and Molecular Data 

All of the capabilities described in the preceding sections 
are only useful to the extent that there are data available 
to describe the level structure, Einstein coefficients, and col- 
lision rate coefficients of various emitter species. To obtain 
these data DESPOTIC relies on the Leiden Atomi c and Molec- 
ular Database (LAMBDA; ISchoier et al.l 120051 ) . which pro- 
vides such data for a variety of species, including the domi- 
nant cooling species in the cold atomic and molecular ISM. 
LAMDA is updated regularly as new calculations or lab- 
oratory measurements of collision rate coefficients become 
available. 

To maximize ease of use, DESPOTIC includes an auto- 
mated downloading and caching capability. When run on a 
machine able to access the internet, DESPOTIC will attempt 
to automatically fetch LAMDA data files as needed. Down- 
loaded files are stored locally for later use. To ensure that 
users take advantage of LAMDA updates as they become 
available, DESPOTIC re-downloads data files that are older 
than six months. A script to force a refresh of the local 
database on a shorter timescale is also included with the 
package. While these capabilities are activated by default, 
users can choose to override them, and manually manage 
their own databases of atomic and molecular data should 
they prefer. 

Finally, DESPOTIC makes three approximations in situ- 
ations where data from LAMDA is not available. First, for 
some species, LAMDA provides estimates only of collision 
rate coefficients for H2, not for 0H2 and pH2 separately, or it 
provides only 0H2 or pHb. In such cases, DESPOTIC assumes 
that the 0H2 and pLb collision rate coefficients are equal, 
and, if only generic H2 rates are given, it sets both of them 
equal to those. 

Second, for some species collision rate coefficients for H2 
are available, but collision rate coefficients for He are not. In 
this case DESPOTIC assumes tha t He collision rate co efficients 
are related to those for H2 by (|Schoier et al.ll2005f ) 



&He = kiln 



Ms-H 2 



1/2 



(66) 



where (j, s -h 2 is the reduced mass of the species s with H2, 
and similarly for /Lt s -He- 

Third, by default DESPOTIC will not extrapolate colli- 
sion rates outside the range of temperatures provided in the 
LAMDA tables. However, the user can override this default 
behavior, in which case DESPOTIC will extrapolate by as- 
suming that the downward collision rate coefficient varies 
as a powerlaw in the gas kinetic temperature. For linear 
molecules, a more accurate extrapolation motivated by a 



io z 



ULIRG 1J CO 
ULIRG 13 CO 
Milky Way 12 C0 
Milky Way 13 CO 




j=l-o .7=2-1 .7=3-2 .7=4-3 .7=5-4 .7=6-5 .7=7-6 .7=8-7 
Transition 



Figure 2. Spectral line energy distribution for the first 8 ro- 
tational transitions of CO and 13 CO, computed for the models 
MilkyWayGMC and ULIRG described in Table[T] The plot shows 
the velocity-integrated brightness temperature in each line nor- 



malized by 7Vh,20 = AW 10 cm ~ 
CMB has been subtracted off. 



The contribution of the 



quantum mechanica l treatment of the collision is possible 
(see the Section 6 of lSchoier et al.ll20050 . but no such treat- 
ment is available for non-linear molecules. Much of what ex- 
trapolation can be done with confidence is already included 
in the LAMDA files. The main reason I include further ex- 
trapolation as an option in DESPOTIC is that it is often neces- 
sary in calculations of equilibrium temperature. This is be- 
cause the automated root-finder, in its attempts to located 
the values of T g and Td for which de gtSp /dt = ded, Bp /dt — 0, 
sometimes explores over a wide range in T g , necessitating 
the calculation of rate coefficients outside the range where 
they are known with confidence. In such cases it is not im- 
portant that the computed rate coefficients be particularly 
accurate, simply that the code not halt when it attempts to 
compute them, and that they not be so wildly inaccurate as 
to cause the root-finder to wander off into unphysical terrain. 
For such purpose, a simple crude powerlaw extrapolation is 
sufficient. 



4 SAMPLE APPLICATIONS 

In this section I provide some sample applications to demon- 
strate DESPOTIC's capabilities. Each of these applications op- 
erates on one or more example clouds, whose properties are 
specified in Table [T] The values given in this Table are in- 
tended to be examples only, but input files corresponding 
to each of them are included with the DESPOTIC library to 
provide example templates that users can modify to set up 
their own clouds. The code to perform each of the example 
calculations listed below is also included with the DESPOTIC 
download. 



4.1 CO Spectral Line Energy Distributions 

As a first example of DESPOTIC's capabilities, Figure[5]shows 
a calculation of CO and 13 CO spectral line energy distribu- 
tions (SLEDs) for the MilkyWayGMC and ULIRG clouds 
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Table 1. Sample clouds 



Cloud Name 




MilkyWayGMC 


ULIRG 




ProtostellarCore 


PostShockSlab 


Physical Propert: 


ies 














n H [cm -3 ] 






10 2 


10 s 




10 2 - 10 8 


10 3 


N H [cm" 2 ] 




l.E 


1 x 10 22 


10 24 




1.0 X 10 23 


1.5 X 10 22 


iTnt [km s -1 ] 






2.0 


80.0 




0.1 


0.5 


T g [K] 






8 


45 




8 


250 


T d [K] 






8 


(50 




8 


8 


Composition 
















ZHI 






0.0 


0.0 




0.0 


0.0 


ZoH2 






0.1 


0.1 




0.1 


0.1 


z p H2 






0.1 


0.4 




0,-1 


0.4 


^He 






0.1 


0.1 




0.1 


0.1 


le 






0.0 


0.0 




0.0 


0.0 


X H+ 






0.0 


0.0 




0.0 


0.0 


Dust Properties 
















CfGD [erg cm 3 K" 


-3/2] 


3.2 


x IO- 34 


3.2 X 10" 


■34 


3.2 X 10" 34 


3.2 X 10~ 34 


°d,iO [cm 2 H" 1 ] 




2.0 


x 10" 26 


2.0 x 10" 


■26 


2.0 x 10~ 26 2.0 x 10~ 26 




°d,PE [cm 2 H" 1 ] 




1.0 


x 10" 21 


1.0 x 10" 


21 


1.0 x 10" 21 


1.0 x 10" 21 


Cd.lSRF [cm 2 IT" 


'] 


3.0 


x 10~ 22 


3.0 x 10" 


22 


3.0 x 10" 22 


3.0 x 10~ 22 


z' d 






1.0 


1.0 




1.0 


1.0 


Pd 






2.0 


2.0 




2.0 


2.0 


Radiation Field Propert: 


ies 












Tcmb [K] 






2.73 


2.73 




2.73 


2.73 


Trad, dust [K] 






0.0 


60.0 




8.0 


8.0 


C [s- 1 H" 1 ] 




1.0 


x 10~ 16 


2.0 x 10" 


-15 


2.0 x 10" 17 


2.0 x 10~ 17 


X 






1.0 


1.0 x 10 4 


1.0 


1.0 


Emitting Species 


Abund 


ances 












CO 




1.0 


x 10" 4 


1.0 x 10' 


-4 


1.0 x 10- 4 


1.0 x 10" 4 


13 C0 




5.0 


x 10" 7 


5.0 x 10" 


-7 


5.0 x 10" 7 


5.0 x 10~ 7 


c 18 o 






- 


- 




5.0 x 10" 8 


5.0 x 10" 8 * 


c 






- 


- 




5.0 x 10" 7 


5.0 x 10" 7 * 


o 






- 


- 




5.0 x 10" 6 


5.0 x 10" 6 


cs 






- 


- 




1.0 x 10" 8 


1.0 x 10" 8 * 


HCO+ 






- 


- 




1.0 x 10" 8 


1.0 x 10" 8 * 


pNH 3 






- 


- 




1.0 x 10" 8 


1.0 x 10" 8 * 


0NH3 






- 


- 




1.0 x 10" 8 


1.0 x 10" 8 * 


pH 2 CO 






- 


- 




1.0 x 10" 8 


1.0 x 10" 8 * 


oH 2 CO 






- 


- 




1.0 x 10" 8 


1.0 x 10" 8 * 


pH 2 






- 


- 




1.0 x 10" 8 


1.0 x 10" 8 * 


oH 2 






- 


- 




1.0 x 10" 8 


1.0 x 10" 8 * 



The table gives initial properties for the example cloud models used in § U For applications where T g and Td are fixed, the values given 

in the table are the values used; for applications where T g and T d are to be calculated, they are used as initial guesses. For the 

ProtostellarCore model, the density is given as a r ange becau s e a ra nge of models are run. The abundances in this model have been 

chosen to roughly match those recommended in Goldsmith (2001). For the PostShockSlab model, emitting species marked with 

asterisks indicate species from which line emission is computed, but that are ignored for the purposes of calculating the thermal 

evol ution. The molecul ar d ata from LAMDA use d in evaluating these mode l s are take n from the fol l owing sour ces: CO, 13 CO, and 

C 18 Q: lYane et all (t2010h ; C: ISchroder et al.l lll99ll) andlStaemmler fc Flowed [|l99lh ; Q :|jaquet et all J1992I); CS:lTurner et alj dl992l) ; 

HCO+: lFlowerl lll999l) ; NH^ lDanbv et all lll988h ; H9CO: lGreenl <199lD ; H,0: |Daniel. Dubernet fc Grosieanl ll201ll) . 



described in Table [T] For this computation, the gas temper- 
ature is left fixed to the input value, and the level popula- 
tions are computing using the escape probability formalism. 
As expected, all lines of the ULIRG are much brighter due to 
its higher gas kinetic temperature and velocity dispersion - 
to first order, the velocity-integrated brightness temperature 
of an optically thick line is simply the product of those two. 
In addition, the falloff in luminosity with J is much slower for 
the ULIRG than for the Milky Way cloud. This is as a result 



of the much higher density and temperature of the ULIRG. 
The former allows its higher levels to be close to thermally 
populated, and the latter causes their thermal populations 
to be large. We also see that the 12 CO(1-0) to 13 CO(1-0) ra- 
tio is larger for the ULIRG than for the Milky Way model, 
reflecting the higher optical depth of the ULIRG. At higher 
J, where the optical depth drops, the line ratios of the two 
isotopomers vary less between the two models. 

Note that this computation for 12 CO (1-0) is equivalent 
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Figure 3. Equilibrium gas and dust temperatures versus density 
for the ProtostellarCore model described in § 14.21 
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to calculating the CO "X-factor" that relates CO intensities 
to cloud masses and column densities. The values calculated 
by DESPOTIC are X C o = 2.6 x 10 20 cm~ 2 / (K km s _1 ) for 
Milky Way CMC, and X C o = 7.0 x 10 19 cm" 2 / (K km s" 1 ). 
This is in line with other theoretical and observational 
estimates for normal galaxies and ULIRGs, respectively 
l|Bolatto. Wolfire fc Lerovll2013h . 

To illustrate the simplicity of the interface, below is the 
entirety of the Python code required to calculate the 12 CO 
and 13 CO SLEDs for the Milky WayGMC model: 

from despotic import cloud 

gmc = cloud(f ileName='cloudf iles/MilkyWayGMC.desp' ) 

gmclines = gmc .lineLum( 'co' ) 

gmclinesl3 = gmc .lineLum( ' 13co' ) 

The first line imports the necessary parts of the DESPOTIC 
library into the local environment. The second reads the in- 
put file that contains the MilkyWayGMC model properties; 
the structure of this file mirrors that of Table [l] The final 
two lines calculate the level populations and line luminosi- 
ties for 12 CO and 13 CO. The code to evaluate the ULIRG 
model is analogous. The total time required to compute all 
four SLEDs (two models, two SLEDs each) is ~ 0.6 seconds, 
not counting the time required to read the input files and 
initialize. 



4.2 Temperatures of Protostellar Cores 

As a second example application, I use DESPOTIC to calcu- 
late the equilibrium gas and dust temperatures in proto- 
stellar cores as a function of density, using the algorithms 
outlines in § 13.2.31 In this calculation I include a large num- 
ber of cooling species (see Table [TJ in order to assess their 
density- and temperature-dependent contribution to cores' 
thermal balance. For this calculation I use the Protostellar- 
Core model in Table [l] I then compute a grid of models 
with densities in the range «h = 10 2 — 10 8 cm -3 in steps 
of 0.2 dex. For each model, I compute the equilibrium gas 
and dust temperatures, and, once the equilibrium has been 
calculated, I record the values of all the heating and cooling 
terms. To illustrate the simplicity of DESPOTIC, note that 
the code required to calculate the equilibrium temperature 
for each model is a single line: 

core . setTempEqO 



Figure 5. Contributions to the overall line cooling rate for in- 
dividual atomic molecular species in the ProtostellarCore model. 
All line cooling rates are computed at the equilibrium temper- 
atures shown in Figure [3] Note that these rates are computed 
for constant abundances, and thus do not properly account for 
depletion at high densities. They are therefore likely to be over- 
estimates at the high-density end, as discussed in § 14.21 



Figure [3] shows the equilibrium temperatures as a func- 
tion of density, Figure [4] shows the contributions of the vari- 
ous heating and cooling processes, and Figure [5] further sub- 
divides the line cooling into the contributions made by indi- 
vidual species. The plots illustrate a number of phenomena. 
First, the gas temperature is relatively high at low densities, 
and drops as the density increases. At densities below ~ 10 4 
cm - this drop is driven by increasingly effective line cool- 
ing. Between 10 4 and 10 5 cm -3 , dust-gas collisions become 
competitive with line cooling, and lock the dust and gas 
temperature together, such that dust-gas energy exchange 
becomes dominant in setting the temperature. The dust in 
turn is always locked close to the infrared radiation field 
temperature, because the IR heating rate and thermal cool- 
ing rate both exceed all other sources and sinks of energy 
for the dust by orders of magnitude. 

In terms of molecular line cooling, at low densities the 
dominant coolants are CO, 13 CO, and C. As the density 
rises and the dust temperature drops, these become less im- 
portant because dust coupling lowers the gas temperature. 
This makes it more difficult to excite the higher J lines 
that have lower optical depths. At the same time, other 
species make an increasing contribution to the cooling as 
the density approaches their critical densities and begins 
to provide efficient collisional excitation. However, I cau- 
tion that these calculations, because they assume constant 
abundances, do not properly model the effects of freeze-out 
onto grain surfaces. This effect will reduce the gas-phase 
abundances of sulfur-bearing molecules at densities above 



10" 
3 



— 10 cm" , carbon-bearing molecules above ~ 10 
cm ", an d nitrogen-bearing mol ecules above ~ 10 7 — 10 s 
cm -3 fe.g. lBergin fc Langerll 19971 ). A more accurate calcula- 
tion could be achieved by using tabulated density-dependent 
abundances. 

This is the most computationally-intensive of the ex- 
ample applications provided, due to the high optical depth 
and the large number of molecular coolants included. The 
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Figure 4. Values of the various heating and cooling terms for dust and gas in the ProtostellarCore models, calculated at the equilibrium 
temperatures shown in Figure [3] The panels show gas heating terms, gas cooling terms, dust heating terms, and dust cooling terms, as 
indicated in the legends. The terms shown are gas ionization heating, r; on , dust IR heating, r^jR, dust ISRF heating r<iiSRF! dust 
line heating T^ i; ne , gas line cooling, Au ne , dust thermal cooling, A,j, and dust-gas energy exchange, Ntgd- The calculations also include 
gravitational and photoelectric heating, but these terms are below the plotted range. 



majority of the computational effort involves iterating to ob- 
tain the level populations at high optical depth. Evaluating 
the entire grid of 31 models requires just over 5 minutes. 
However, since only a few chemical species are actually im- 
portant to the thermal balance, one could obtain the results 
far more quickly simply by ignoring the large number of 
energetically-unimportant species when calculating the tem- 
perature, and only calculating their line luminosities once 
the temperatures have converged. DESPOTIC includes a capa- 
bility to mark certain species as energetically-unimportant, 
allowing them to be treated in precisely this manner, and I 
demonstrate this capability in the next example. 



4.3 Time-Dependent Cooling of Post-Shock Gas 

A third example, which makes use of DESPDTIC's ability to 
calculate time-dependent temperature evolution (£j |3.2.3[) . is 
to calculate the cooling of out-of-equilibrium gas. I consider 
a slab of gas whose properties are given by the PostShock- 
Slab model in Table [T] At time t — 0, the gas has just 
been shock-heated to an out-of-equilibrium temperature of 
250 K, and I calculate the time evolution of its tempera- 
ture and line emission thereafter, assuming that the gas is 
isobaric and using a slab geometry to compute escape prob- 
abilities. In calculating the thermal evolution I include only 
the energetically-dominant coolants CO, CO, and O, but 
I also periodically compute the line emission of a large num- 
ber of other species as well. By making this assumption, the 
total computer time required to evolve the model 40 kyr, 
including periodic calculation of emission from many lines, 
is ~ 10 minutes. 

Figure [6] shows the gas temperature, dust temperature, 
and gas density versus time as computed by DESPOTIC for 
this initial condition. Figure shows the contributions of 
various species to the cooling. As the plot shows, cooling 
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Figure 6. Gas temperature, dust temperature, and gas density 
versus time for isobaric cooling of the PostShockSlab model, as 
described in § 14.31 



is dominated by CO lines, with minor contributions from 
13 CO, O, and dust, and negligible contributions from all 
other sources. In Figure [5] I further examine the cooling, 
by showing how the CO spectral line energy distribution 
changes with time. As the plot shows, the SLED initially 
peaks near J = 7 — 6, and moves to a cooler SLED at time 
passes. At the final time shown, J — 3 — 2 is the dominant 
coolant. Note that this differs from the result shown in Fig- 
ure [2] for a typical CMC because the post-shock slab we are 
considering has a significantly lower velocity dispersion and 
a significantly higher density. Both of these favor cooling 
through higher J lines, the former because it increases the 
optical depth for low J lines, and the latter because it helps 
to thermalize higher J states. 
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Figure 7. Rates of cooling provided by CO lines, 13 CO lines, O 
lines, and dust versus time, for the PostShockSlab model shown 
in Figure [6] The gray thick line shows the sum of all coolants, 
including all a number of lines that are not shown because they 
all below the range of cooling rates plotted. 



Figure 9. Brightness temperature versus velocity relative to 
line center for the lines CS(3 — 2) and C 34 S(3 — 2) produced by 
a collapsing protostellar core. The contribution of the CMB has 
been subtracted off. Details of the core parameters are given in 
§111 
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Figure 8. CO line spectral line energy distribution for the Post- 
ShockSlab model shown in Figure |S] Each of the lines shows the 
relative contributions of the indicated rotational transitions of 
CO to the total cooling rate at the indicated times of kyr, 20 
kyr, and 40 kyr. Contributions are normalized so that the sum 
over all transitions is unity. 



4.4 Inverse P Cygni Profiles 

As a final application, I use DESPOTIC to calculate line pro- 
files in a collapsing protostellar core. For this example, I 
consider a core with a radius of R — 0.02 pc with a velocity 
profile v(r) = — 0A(r/R)f km s _1 . The temperature profile 
is T(r) = 8 + I2exp(-2r 2 /7? 2 ) K, so that the temperature 
reaches a peak of 20 K at the center, dropping close to 8 
K at large radii. The core also has a position-independent 
non-thermal velocity dispersion of 0.2 km s . 

For this core I use DESPOTIC to compute the pro- 
files of the CS(3-2) and C 34 S(3-2) lines. The former is 
frequently used to measur e inverse P Cygni profiles (e.g. 
iLee. Myers fc Plume! 120041 ). while the latter is its low- 
abundance isotopomer, and thus should remain optically 
thin and not show inverse P Cygni profiles, f take the den- 
sity of CS molecules in the core to be 0.1 cm -3 , independent 
of position, and the abundance of C 34 S to be 1/22 of this 



value, corresponding to the abundance ratio of S to S 
measured on Earth. 

Figure [9] show the line profiles computed by DESPOTIC. 
As expected, the marginally optically thick CS line produces 
a double-peaked asymmetric inverse P Cygni profile, indica- 
tive of infall. The lower abundance C 34 S is optically thin 
and produces a symmetric profile of lower total intensity. 
The total time required to perform the computation is ~ 10 



5 LIMITATIONS AND CAVEATS 

While DESPOTIC provides reasonable estimates of the ther- 
mal behavior and spectra of interstellar clouds over a wide 
range of environments, it also has significant limitations, 
which 1 discuss here as a warning to potential users. The 
major limitations of the code are: 

• DESPOTIC's treatment of dust temperatures is very 
crude in the regime of clouds that are optically thick to 
their own cooling radiation. In such clouds the dust tem- 
perature will be determined largely by the value of Trad, dust 
that the user selects, ff a user requires accurate dust tem- 
peratures i n such clouds, he or sh e is advised to use a code 
like dusty (llvezic fc Elitzurlll997h to calculate the dust tem- 
perature and radiation field within the cloud, then use this 
to set T^d. dust for the purposes of a DESPOTIC calculation. 

• DESPOTIC neglects the contribution of the dust radia- 
tion field to the photon occupation number when calculating 
level populations, on the grounds that, because dust opti- 
cal depths are small at low frequencies, such fields are of- 
ten highly sub-thermal at the low frequencies where most 
important molecular lines lie. However, in some circum- 
stanc es, e.g. protostellar disks IjKrumholz. Klein fc McKeel 
120071 ). the column density is so high that dust optical depths 
can exceed unity even at frequencies as low as 20 GHz. In 
such environments excitation and de-excitation of molecules 
by interaction with the infrared field is non-negligible, and 
DESPOTIC will not give accurate results. 

• DESPOTIC uses a one-zone model, and this is not capa- 
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ble of capturing effects that depend on radiative transfer. 
In particular, DESPOTIC cannot handle maser emission, and 
it cannot handle effects on the line shape that arise from 
spatially-variable departures from LTE. 

• DESPOTIC does not perform any chemical calculations. 
It is up to the user to input the chemical abundances, and 
the results DESPOTIC produces will only be as good as those 
abundances. More subtly, it does include the effects of selec- 
tive chemical destruction of excited states on line emission, 
and it does not include any heating or cooling of the gas 
or dust as a result of chemical reactions, such as heating 
of dust gr ains by exothermic formation of H2 on grain sur- 
faces (e.g. lLesaffre et al. 2005). DESPOTIC provides a mecha- 
nism to include chemical heating and cooling, since the user 
can specify arbitrary additional heating and cooling terms, 
but it is up to the user to determine whether there are any 
energetically-important chemical reactions for the problem 
under consideration, and, if so, to implement the necessary 
code. 



6 SUMMARY 

I introduce DESPOTIC, a Python-based, open-source software 
library for calculating spectra, heating and cooling rates, 
and time-dependent and time-independent thermal proper- 
ties of optically thick interstellar clouds. DESPOTIC includes 
all the dominant heating and cooling processes for both gas 
and dust over a wide range of interstellar environments, and 
can be used to conduct both fast sweeps of parameter space 
and interactive explorations within an interactive Python 
environment. It is intended to allow theoretical investiga- 
tors to obtain approximate values of parameters such as 
cloud temperatures, major heating and cooling processes, 
and observable line emission, without the difficulty and time 
investment of developing their own statistical and thermal 
equilibrium codes, and with significantly less investment of 
CPU and human time than would be required to approach 
such problems using a detailed PDR code. DESPOTIC is under 
continued development, and additional features capabilities 
will be released to the community as they are implemented. 
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APPENDIX A: SPECIFIC HEATS 

Calculating the time evolution of the temperature requires 
knowing the specific heat per H nucleus at constant volume 
c„,h, defined by 



£Vh = 



J_(deg\ 



(Al) 



n H KdTjp' 

where e g is the gas internal energy per unit volume, given 
by 

, d In z a 



e g — 2 y risksT 



dlnT' 



(A2) 



where the sum runs over species s, n s is the number density 
of species s, and z s is the partition function per unit volume 
for that species. The latter is given by 

Zs — ^s. trans^s, rot ^s,vib-^s, spin , V / 

where the terms appearing in the equation above are the par- 
tition functions for the translation, rotational, vibrational, 
and spin degrees of freedom of species i. In principle we 
should also include a term describing electronic degrees of 
freedom, but at the relatively low temperatures for which 
DESPOTIC is intended, we can safely assume that these are 
not excited. For all the species included in DESPOTIC except 
molecular hydrogen (i.e. for H I, He, H + , and e), the con- 
tribution of the specific heat is trivial, because all of the 
partition functions except translation and spin are unity, 
and the spin term is temperature-independent. Thus for all 
these species 

dlnz s _ <91n/7 Sitr an S _ 3 
dhiT ~ <91nT ~~ 2' 



(A4) 



For ortho- and para-H2 on the other hand, Z TO t and Z v jb are 
not unity dBlack fc Bodenheime"r] 1 19751 ; iBolev et all 120071 ; 
iTomida et al.ll2013n : 



•^oH 2 ,rot 

•^pH 2 ,rot = 
■^H2,vib = 


= 2^ 3 ( 2J + 1 ) ex P 

J odd 

= ^2 ( 2J + !) ex P 

J even 

l 


1 - exp(-6» vib /T) 



J(J+l)6> r , 
T 

J {J + l)0 ro i 



T 



(A5) 

(A6) 
(A7) 



where 6 TOt = 85.3 K and 8 vih = 5984 K. Note that the vi- 
brational partition function is the same for ortho- and para- 
H2, but the rotational partition functions are different. With 
these partition functions, the energy per unit volume includ- 
ing all species is 



kn 



-T 2_^n s +n p n 2 



pH'2 ! r °t 



\ZpH 2 



t)T 



+ n n 2 



( T 2 8Z oK2 



<:)T 



+ (n p H 2 + f»oH 2 ) Ovib 



exp(-g vib /T) 
1 - exp(-0 vib /T) ' 



(A8) 



where again the sum runs over all all species. 

Deriving the specific heat c v from this expression re- 
quires making an assumption about how the number den- 
sities of ortho- and para-H2 vary with temperature. At the 
low temperatures found in interstellar clouds, there is gener- 
ally no efficient mechanism for converting between the two 
states, and thus the most reasonable assumption is that 
these number densities are temperature-independent. Ob- 
servations showing that the ortho- to para- ratio in molec- 
ular clouds is far from equilibrium (e.g. iNeufeld et alj|2006l : 
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iPagani. Roueff fc Lesaffrel l201ll ; iDislaire et all I2012T ) sup- 
port this assumption. For temperature-independent values 
of n P H 2 and n H 2 , we therefore have 



Cy.H 

k B 



3 r , d ( T 2 0Z pil2 



+ X Uo 



dT 



T 2 dz oW . 

?oh, dT 



+ (x p H 2 + ZoHa 



6>; ib exp(-e vib /r) 

T 2 [l-exp(-0 vib /T)p 



(A9) 



Note that this expression involves the abundances ratios x 
rather than number densities n because we have normalized 
all quantities to the number density of H nuclei. The specific 
heat at constant pressure is simply c Pi H/fcs = c V: u/kB + 1. 



APPENDIX B: COMPARISON BETWEEN 

DESPOTIC AND RADEX 

As a check on DESPOTIC and to illustrate its strengths and 
weaknesses, in this Appendix I pr ovide a detailed compar i- 
son between DESPOTIC and RADEX (jvan der Tak et al.ll2007h . 
RADEX does not include DESPOTIC's capabilities for com- 
puting heating and cooling rates, thermal equilibria, time- 
dependent thermal evolution, or line shapes, so this test is 
limited to the capabilities that the two codes have in com- 
mon: computing level populations and emergent line inten- 
sities from a cloud of specified physical properties. 

For the purposes of this test, I compute the CO spec- 
tral line energy distribution for a cloud with temperature 
of T g — 10 K, a full-width-at-half-maximum velocity spread 
of 2.0 km s _1 , and a CO abundance :rco = 10~ 4 over a 



grid of volume densities from nn 



-2 



10 cm and col- 



umn densities Na = 10 14 — 10 /4 cm~' ! , in steps of 0.2 dex in 
both dimensions. The grid is chosen to cover a wide range of 
conditions, from optically thin to optically thick, and from 
thermalized (for the first for levels) to highly sub-thermal. 
For both codes the background radiation field is set to the 
CMB value of 2.73 K, and I use slab geometry for the es- 
cape probability calculation, since RADEX and DESPOTIC use 
the same approximate expression for the escape probabil- 
ity in that case. I perform the RADEX computation using a 
slightly modified version of the radex_grid . py wrapper that 
is distributed as part of the RADEX package. 

To ensure that the computations are identical, for the 
DESPOTIC calculation I set the non-thermal velocity disper- 
sion to a temperature-dependent value ont in DESPOTIC to 
ctnt = [FWHM 2 /81n2 + c^/^com H ] 1/2 , where /xco = 28 
is the molecular weight of CO. This guarantees that the 
velocity dispersions are the same in the two calculations. 
Similarly, I disable clumping and I set all dust opacities to 
in DESPOTIC, since RADEX includes neither clumping nor 
dust absorption. Finally, I set the abundances of ortho- and 
para-H 2 in DESPOTIC to 



2-oHo — 



9e 



-2e rot /T 9 



1 _|_ 9 e -26>rot/T g 
1 



(Bl) 
(B2) 



pH2 ~~' i + gc-aWi' 
consistent with RADEX's hardwired assumption that the ratio 



of ortho- to para-H2 is given by the thermal ratio of the 
populations of the H2 J = 1 to J = states. 

Comparison of the results indicates that the level pop- 
ulations and line optical depths returned by the two codes 
are identical to the level of precision with which RADEX 
writes output. The line fluxes returned by the codes, inter- 
estingly enough, are not identical, and this is due to a minor 
lack of self-consistency in the escape probability approxima- 
tion itself. To compute the frequency-integrated line flux, 
DESPOTIC first computes the total rate of energy emission 
per H nucleus from equation (|41|) . and then computes the 
integrated intensity and brightness temperature from equa- 
tions (|60p and (|61[) . In contrast, RADEX computes the output 
intensity using the transfer equation for a uniform medium. 
It uses the level populations to compute an excitation tem- 
perature T GXj ij between every pair of levels i, j, computes 
the optical dep th at line center rg from the level populations 
(equation 21 of Ivan der Tak et al.ll2007l ). and then computes 
the emergent integrated intensity as 

/ [B v {Tk,h) (l - e~ T ^) + e~ T ^ B„(T C mb)] 4> v dv, (B3) 

where (j> v is the line shape function, which is taken to be 
a Gaussian whose dispersion is determined by the input 
FWHM. To obtain the cooling rate per H nucleus, this quan- 
tity is simply divided by the total column density. In the 
limit of high optical depth, and neglecting the contribution 
of the background radiation field (which is indeed negligible 
in the example given) , with some algebra one may show that 
RADEX's expression reduces to 



A., 



-AijAEijfi, 



(B4) 



while DESPOTIC's expression (equation 14 1 p reduces to 

A B ,ij = Aj AijAEijfi. (B5) 

In the optically thin limit, the factors of 1/Ty and /3ij are 
omitted, rendering the expressions identical. As expected, in 
the optically thin limit the two codes produce results that 
are identical to the precision with which RADEX writes out- 
put. For optically thick lines, on the other hand, the two ex- 
pressions above are identical only if fiij — >• 1/tjj as r»j — > 00. 
This is the case for the LVG approximation, and for the ex- 
pression RADEX uses for spherical geometr^Jj, but it is not true 
for slab geometry or for the approximation that DESPOTIC 
uses in spherical geometry. 

This disagreement arises from a fundamental limitation 
of the escape probability approximation. In this approxima- 
tion, one assumes that there is a uniform escape probabil- 
ity that characterizes the entire cloud, and that the level 
populations within the cloud are also uniform, but these 
two assumptions are not fully consistent. DESPOTIC's cal- 
culation of the line luminosity follows from the former as- 
sumption, while RADEX's follows from the latter. However, 
the former assumption is preferable for the types of prob- 
lems that DESPOTIC is intended to solve, because it enforces 



9 In spherical geometry one must be careful to correct for the 
fact that RADEX defines the optical depth appearing in equation 
l!B4t as measured along a cloud diameter, which is larger than the 
projection-averaged optical depth used in DESPOTIC by a factor of 

3/2. 
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strong consistency between the rate of photon emission and 
escape from the cloud, and the rate of rate of energy loss via 
line cooling. 

Finally, I note that, in timing tests, RADEX performs 
this calculation a factor of ~ 5 faster than DESPOTIC. This 
difference is not surprising, given that RADEX is a single- 
purpose tool written in Fortran, compiled with heavy op- 
timization, and where many decisions are made at compile 
time (e.g. the geometry used to compute escape probabili- 
ties), while DESPOTIC is a much more general-purpose and 
interactive tool written in a non-optimized language, and 
with a large number of options that are specified at run- 
time. In neither case is the computational cost prohibitive, 
however. On the workstation where I performed the tests, 
the full grid of 1581 models required roughly 35 seconds to 
evaluate for RADEX, and a bit under 3 minutes for DESPOTIC. 



APPENDIX C: HANDLING 
ILL-CONDITIONED MATRICES 

As discussed in the main text, for some species the ma- 
trix M that appears in equation (|57p can be extremely ill- 
conditioned. Figure [Cll illustrates the nature of the problem 
by graphically displaying M for the species CO, C + o and 
0NH3, all computed for a cloud with nu = 10 3 cm -3 com- 
posed of pure pH2 at a temperature of 10 K, embedded in the 
cosmic microwave background at temperature Tcmb = 2.73 
K. For simplicity the cloud is assumed to be optically thin 
and to have a clumping factor f c — 1. The matrix describing 
CO has a wide range of transition rates, but every row and 
column contains at least one transition rate coefficient whose 
magnitude is comparable to that of the largest elements of 
the matrix. As a result, the CO matrix is well-conditioned: 
for the example shown in Figure [Cll the condition number is 
135. This presents no challenges for numerical solution. On 
the other hand, the matrices for C + and 0NH3 both have 
the property that a few elements are much larger than most 
of the rest of the matrix. As a result they have condition 
numbers of 2.7 x 10 36 and 1.2 x 10 13 , respectively. These 
condition numbers imply that a numerical solution to equa- 
tion (157[1 for C + and 0NH3 would have ~ 36 and ~ 13 fewer 
digits of accuracy than machine accuracy, rendering numer- 
ical solutions obtained for these matrices meaningless. 

The high condition numbers of the matrices are a di- 
rect result of the physical processes they describe, and the 
divergence in timescales between transitions between differ- 
ent states. In the matrix for C + , for example, the largest 
elements correspond to transitions between the 2 Si/2 and 
2 D 3 / 2 states (states 7 and 6, respectively, in Figure [Cl) and 
the ground, 2 Pi/ 2 state (state in Figure [CT|) . These have 
Einstein coefficients A ~ 10 9 s _1 , compared with the fine- 



structure transition between the first two states, which has 
A = 2.3 x 10 -6 s _1 . Similarly, for 0NH3, the largest elements 
of M describe transitions such as (J, K) v = (7, 6)1 — > (6, 6)0 
and (7, 6)0 — > (6, 6)1 (elements ij — 11, 18 and 12, 17, respec- 
tively, in Figure [Cll) . with A ~ 0.1 s~ , while the inversion 

10 LAMDA offers two data tables for C + , one including only the 
low-lying fine-structure levels, and one also including the higher 
energy levels connected to them by UV lines. For the purposes of 
this example, I use the data file including the UV levels. 




Figure CI. A graphical representation of the matrices M cal- 
culated for CO, C+, and 0NH3 using the conditions described 
in Appendix \C\ The color of each block represents the value of 
the corresponding element ij of M, excluding the elements with 
i = N + 1, which are all unity. The color scale is normalized 
and logarithmic, so that the largest element of M is shown in 
red, while dark blue corresponds to a value of 10 — 15 times the 
value of the largest element. Values on the diagonal are masked, 
since they are negative. Recall that the value of an element Mij 
is the sum of the rate coefficients describing transitions into state 
i from state j (including both collisional and radiative transi- 
tions), normalized by the sum of the rate coefficients out of state 
i into any other state. Thus elements above the diagonal represent 
downward transitions, while those below the diagonal are upward 
transitions. 



transitions that are most commonly observed (e.g. (6,6)1 — > 
(6,6)0, element 11,12), have A ~ 10~ 7 s _1 . 

DESPOTIC implements two strategies when it encounters 
an ill-conditioned matrix. First, in many cases high condi- 
tion numbers are associated with large rates for downward 
transitions from high-energy levels. In the simple one-zone 
statistical equilibrium model used by DESPOTIC, the popu- 
lation of any level will be bounded between the values ex- 
pected when the atom is in LTE at T g and when it is in 
LTE at Tcmb- DESPOTIC calculates these two limiting val- 
ues, and if it finds that they are below a numerical flooj 11 !. 
it simply sets the populations of those levels to the floor, 
and removes the associated rows and columns from matrix 
M. If these rows and columns contain large elements, the 
condition number of the matrix is likely to be reduced. For 
the examples shown in Figure IC1I applying this procedure 
eliminates the 6 highest energy levels for C + (and thus the 
six bottom- and right-most rows and columns in M) and the 
11 highest energy levels for 0NH3. In turn, this reduces the 
condition number for the C + matrix to 7.5 x 10 4 , low enough 
to allow numerical solution with tolerable accuracy. Figure 
C2I shows the same graphical representation of the matrix 
for C + as in Figure [Cll after this level reduction procedure. 

Unfortunately this procedure alone still leaves the ma- 
trix for 0NH3 with an unacceptably-high condition number 
of 7.6 x 10 11 . This is because not all of the large matrix ele- 
ments for this case apply only to the high-energy levels. The 
levels (3,3)o, (3,3)i, (4,3)o, and (4, 4)0 are close enough to 
the ground state in energy for their populations not to be 
entirely negligible, but they still have Einstein A values far 
larger than those associated with the inversion transitions. 

The second strategy DESPOTIC employs is to eliminate 
levels whose populations will be small because the rate coef- 
ficients for transitions out of them greatly exceeds the rate 



11 DESPOTIC sets this floor equal to the machine epsilon value for 
the platform on which it is operating, which is usually ~ 10 — 15 . 
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oNH, 




Figure C2. Same as Figure fCll for C + and 0NH3, after level 
reduction as described in the text. 



coefficients for transitions into them. Specifically, consulting 
equation (|35[1 . the total rate at which a particle in state i 
will transition to another state is 



t i.ou 



E 

k 



qik + Pik (1 + n-y : ik) Aik + Pki — n^^iAki 

9i 



.(CI) 



The total rate of transitions into state i is given by the left- 
hand side of equation (|35[1 , and since each /_, is strictly less 
than unity, the rate of transitions into state i is bounded 
above by 



r -<E 



9j 



(C2) 



Thus the population of state i is strictly bounded above by 
r • 



Ji < fi,lm 



-L i.out 



(C3) 



For ill-conditioned matrices, this ratio is very small for some 
levels and very large for others. In the example of 0NH3, once 
the high-temperature levels have been eliminated, the value 
of /i.iim runs from a minimum of 1.3 x 10~ to a maximum of 
1.2 x 10 5 . Thus if M remains ill-conditioned after the high- 
temperature levels have been eliminated, DESPOTIC finds the 
level with the smallest value of /i,ii m and eliminates it in 
exactly the same manner as the high temperature levels. If 
necessary it repeats this procedure with the next smallest 
value of /ijim and so forth, until the condition number of 
the matrix is acceptably small. Figure [C2l shows the matrix 
for 0NH3 after this procedure is complete, leading to the 
elimination of all levels by the five lowest energy ones. The 
condition number of the resulting matrix is 3.5 x 10 5 , again 
allowing accurate numerical solution. 
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